/*
 *  lib-smv: Matrix - Vector Library targeted to Structural Mechanics Problems
 *  Copyright (C) 2006-2008  St.Brcic, Lj.Zugic
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
/** @file Arrays.h
 *  @brief Definition of Vector and Matrix classes of double numbers
 */

#ifndef _ARRAYS_H
#define _ARRAYS_H

#include <iostream>
#include <fstream>
#include <iomanip>
#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <cassert>
#include <string>
#include "iArrays.h"

using namespace std;

namespace smv_lib
{

//! default size of vector/matrix
const int M_DEFAULT = 2;
//! small number - equal to zero (for LU decompositions)
const double TINY = 1.0e-20;  
//! small number - equal to zero (for printing vectors and matrices)
const double SMALL = 1.0e-15;  

//---------- forward declaration
class Matrix;

//-------------------------------------------------------------------------
/**
** \brief class Vector of double numbers
**
**  VECTOR: VECTOR class of double numbers
**  Constructors: 
**  (1) empty (default) constructor: (size is assigned as default value 2, 
**      name is given as "VECTOR", but data should be given)
**  (2) constructor with a size 'n' and all elements 'd' (name is not given)
**  (3) constructor with a size 'n' and a given Name 
**  (4) constructor form a vector file (5) copy constructor
**/
//-------------------------------------------------------------------------
class Vector
{
  private:
    //! good for Matrix-Vector manipulation
    friend class Matrix;
    //! array of double values
    double *data;
    //! number of vector elements
    int len;
    //! input stream
    ifstream in;
    //! output stream
    ofstream out;
    //! name of the Vector
    string Name;
    //! true if the Name is given
    bool isNameGiven;

  public:
    //! empty (default) constructor
    Vector();
    //! Vector with 'n' double elements of value 'd'
    Vector( int n, double d = 0.0 );
    //! Vector with given length 'n' and a name (all elements set to zero)
    Vector( int n, string name );
    //! construct a Vector from a file
    Vector( const char* file );
    //! copy constructor
    Vector( const Vector& original );
    //! copy constructor: convert from int iVector
    Vector( const iVector& original );
    //! destructor
    ~Vector();

    //! assignment operator
    Vector& operator = ( const Vector& original );
    // -------------------- get and set vector elements
    //! set an element (zero based, operator form)
    double& operator[] ( int i );
    //! get an element (operator form)
    const double& operator[] ( int i ) const;
    // -------------------- get and set functions (alternative)
    //! get number of elements
    const int getLen() const { return len; }
    //! get element no 'i' (zero based)
    const double getData( int i ) const;
    //! get the name
    const string getName() const { return Name; }
    //! set length of a Vector
    void setSize( int n ) { len = n; }
    //! set element no 'i' with double data 'd'
    void setData( int i, double d );
    //! set the name
    void setName( string name );

    //------------------- basic linear algebra - operator functions 
    //                    (better use equivalent functions)
    //! Operator '+' addition with a vector: b = self + v
    Vector operator + ( const Vector& v );
    //! Operator '-' subtraction of a vector: b = self - v
    Vector operator - ( const Vector& v );
    //! Operator '*' product with a scalar: v = self * s;
    Vector operator * ( double s );
    //! Operator '*' scalar product with a vector: d = self_Trans * v
    const double operator * ( const Vector& v ) const;
    //! Operator '+=' addition with a vector: self = self + v
    Vector& operator += ( const Vector& v );
    //! Operator '-=' subtraction of a vector: self = self - v
    Vector& operator -= ( const Vector& v );
    //! Operator '*=' product with a scalar: self = self * s
    Vector& operator *= ( double s );

    //! is self == v ? (true if equal, the names are not important)
    bool operator == ( const Vector& v );

    //------------------- basic linear algebra - equivalent functions
    //! addition with a vector: b = self + v
    Vector Add( const Vector& v );
    //! subtraction of a vector: b = self - v
    Vector Sub( const Vector& v );
    //! product of the vector with scalar:  v = self * s
    Vector Dot( double s );
    //! dot product with another vector: d = self_Trans * v
    const double Dot( const Vector& v ) const;
    //! outer product with a vector: B = self * v_Trans 
    Matrix OuterDot( const Vector& v );
    //-------------------- subVector manipulations
    //! Add sub-vector into the vector (add to the old content)
    Vector& addSubVector( int i, const Vector& Sub );
    //! Take sub-vector from a given vector
    const Vector getSubVector( int i, int subLen );

    // ------------------- arithmetic operations
    //! change the size of a Vector
    void Resize( int n );
    //! make all elements equal to zeroes
    void Zero( void );
    //! max element in a Vector
    const double Max( void ) const;
    //! min element in a Vector
    const double Min( void ) const;
    //! mean value of all elements
    const double Mean( void ) const;
    //! clean vector of small numbers
    void Clean( void );
    //-------------------- norms of the vector
    //! norm of Self vector (intensity): (Self * Self)^(1/2)
    const double Norm( void ) const;
    //! B_Norm of Self vector: (Self_T * B * Self)^(1/2)
    const double B_Norm( const Matrix& B ) const;

    // ------------------- console output
    //! print all in one row
    void rPrint();
    //! print each element in a separete line
    void Print();
    //! ------------------- file I/O
    //! read from predefined file "Name.txt"
    Vector Read();
    //! write to predefined file  "Name.txt"
    void Write();
    //! read from a filev
    Vector ReadFrom( const char* file );
    //! write to a file (L entries in a line)
    void WriteTo( const char* file, int L = 20 );
    // ------------------- stream I/O
    //! output stream
    friend ostream& operator << ( ostream& s, const Vector& v );
};

//================================== EXTERNAL VECTOR FUNCTIONS
//
//! addition of two vectors: v3 = V1 + v2
Vector operator + ( const Vector& v1, const Vector& v2 );
//! subtraction of two vectors: v3 = v1 - v2
Vector operator - ( const Vector& v1, const Vector& v2 );
//! scalar (dot) product of two vectors: d =v1 * v2
double operator * ( const Vector& v1, const Vector& v2 );
//! scalar (dot) product of two vectors: d = v1 * v2 (equivalent form)
double Dot( const Vector& v1, const Vector& v2 );
//! outer product of two vectors: the result is matrix B = v1 * v2_Trans
Matrix OuterDot( const Vector& v1, const Vector& v2 );
//! product of a vector and a double scalar:  v1 = x * v
Vector operator* ( double x, const Vector& v );
//! product of a vector and a double scalar:  v1 = v * x
Vector operator* ( const Vector& v, double x );
//! Euclid norm of a vector v: (intensity of the vector)
double Norm( const Vector& v );
//! B_norm of a Vector: Matrix norm (v_T * B * v)^(1/2)
double B_Norm( const Matrix& B, const Vector& v );
//! true if two vectors are equal (the names are not important)
bool areEqual( const Vector& v1, const Vector& v2 );
//! output stream for a vector
std::ostream& operator << ( ostream& s, const Vector& v );

//----------------------------------------- end of Vector class

//================================= MATRIX CLASS - ONLY DENSE MATRICES

//--- forward declaration of Crout and Cholesky classes
class Crout;
class Cholesky;

//! useful enumerations for matrix classification
enum MAT_TYPE {NOT_CLASSIFIED, DIAGONAL, TRIDIAG, TRI_SYMM, SYMMETRIC, GEN_TYPE};
// diagonal, tridiagonal, tridiagonal + symmetric, symmetric and general
//

//-----------------------------------------------------------------------
/**
** \brief class Matrix of double numbers
**
**    MATRIX: MATRIX class of double numbers
**
**    Only dense matrices are considered
**    Special structure of a matrix (if any) is taken into account:
**    matrices may be diagonal, tridiagonal, symmetric or general
**/
//-----------------------------------------------------------------------
class Matrix
{
  private:
    //! Matrix data (Vector of rows m)
    Vector* m;
    //! Number of rows and columns
    int rows, cols;
    //! Input stream
    ifstream in;
    //! Output stream
    ofstream out;
    //! Name of the Matrix
    string Name;
    //! Indicator: true if the Name is given, false for default name "MATRIX"
    bool isNameGiven;
    //! Matrix type (diagonal, tridiagonal, symmetric or general)
    MAT_TYPE Mat_Type;
    //! Indicator: true is matrix is classified
    bool isClassified;
    //! Indicator: true is matrix is singular (false for regular matrix)
    bool isSingular;
    //! Indicators: diagonal, tridiagonal, symmetric or general
    bool isDiagonal, isTridiagonal, isTriSymmetric, isSymmetric, isGeneral;

    //! Matrix inversions and solutions of linear equations
    //! Inverse of the matrix A
    Matrix* A_inv;
    //! General matrix - Crout's LU decomposition
    Crout* m_Crout;
    //! Symmetrix matrix - Cholesky decomposition
    Cholesky* m_Chol;
    //! Indicator: true if Crout LU is selected
    bool isCrout;
    //! Indicator: true if Cholesky LLT is selected
    bool isCholesky;
    //! Determinant of the matrix A
    double Det;

  public:
    //----------------------- constructors
    //! empty (default) constructor
    Matrix();
    //! Matrix with given rows, columns and all elements set to a double 'd'
    Matrix( int r, int c, double d = 0.0 );
    //! Matrix with given rows, columns and a name (all elements set to zero)
    Matrix( int r, int c, string name );
    //! construct a Matrix from a file
    Matrix( const char* file );
    //! copy constructor
    Matrix( const Matrix& M );
    //! copy constructor: convert from int iMatrix
    Matrix( const iMatrix& M );
    //! destructor
    ~Matrix();
    //! assignment operator
    Matrix& operator = ( const Matrix& M );

    //---------------------- get functions
    //! get number of rows
    int getRows() const { return rows; }
    //! get number of columns
    int getCols() const { return cols; }
    //! get element (i,j) - zero based
    const double getElement( int i, int j ) const;
    //! get elements in a row 'i' of a Matrix
    const Vector getRowData( int i ) const;
    //! get elements in a column 'j' of a Matrix
    const Vector getColData( int j ) const;
    //! get the inverse matrix
    const Matrix getInv();
    //! get determinant of the matrix
    const double getDeterminant();
    //! Indicator: true if inverse matrix is determined
    bool isInverse;
    //! get the name of a Matrix
    const string getName() const { return Name; }
    //! get the matrix type
    MAT_TYPE get_MatType() const { return Mat_Type; }
    //! give written info about the matrix type
    string give_MatType() const;

    //--------------------- set functions
    //! set number of rows
    void setRows( int r ) { rows = r; }
    //! set numer of columns
    void setCols( int c ) { cols = c; }
    //! set the value 'val' to row 'i' and column 'j'
    void setElement( int i, int j, const double val );
    //! set all data in a row 'i' with Vector 'v'
    void setRowData( int i, const Vector& v );
    //! set all data in a column 'j' with Vector 'v'
    void setColData( int j, const Vector& v );
    //! set a name to the Matrix
    void setName( const string name );
    //! set the matrix type
    void set_MatType( MAT_TYPE type ) { Mat_Type = type; }
    //! set matrix classification indicator
    void setIsClassified(bool set) {isClassified = set;}

    //---- basic matrix manipulations
    //! the sum of elements in a row
    const double RowSum( int i );
    //! the sum of elements in a column
    const double ColSum( int j );
    //! matrix is assigned with all zero elements
    void Zero( void );
    //! matrix becomes a unitary matrix
    void Identity( void );
    //! matrix becomes diagonal (all eleemnts are the same: 'd')
    void Diagonal( double d );
    //! matrix becomes diagonal (elements are equal to Vector 'v')
    void Diagonal( const Vector& v );
    //! maximum element on main diagonal
    const double MaxDiagElem( void );
    //! maximum element in whole matrix (square or rectangular)
    const double MaxElement( void );
    //! clean matrix of small numbers
    void Clean( void );
    //------------------------- tridiagonal matrix
    //! true if tri-diagonal matrix is transformed into three vectors a,b,c (b is diagonal)
    bool fromTriDiag( Vector& a, Vector& b, Vector& c );
    //! Three vectors a,b,c are transformed into tri-diagonal matrix A
    Matrix setTriDiag( Vector& a, Vector& b, Vector& c );
    //------------------------ subMatrix manipulations
    //! Insert sub-matrix into the matrix (replace the old content)
    Matrix& insertSubMatrix( int i, int j, const Matrix& Sub );
    //! Add sub-matrix into the matrix (add to the old content)
    Matrix& addSubMatrix( int i, int j, const Matrix& Sub );
    //! Take sub-matrix from a given matrix
    Matrix getSubMatrix( int i, int j, int subRow, int subCol );

    //---- arithmetic operators: SELF = A
    // addition and subtraction
    //! Addition: Matrix C = Matrix A + Matrix B
    Matrix operator + ( const Matrix& B );
    //! Subtraction: Matrix C = Matrix A - Matrix B
    Matrix operator - ( const Matrix& b );
    //! Addition: Matrix A = Matrix A + Matrix B
    Matrix& operator += ( const Matrix& B );
    //! Subtraction: Matrix A = Matrix A - Matrix B
    Matrix& operator -= ( const Matrix& b );
    // various products
    //! Product: Matrix B = Matrix A * double s
    Matrix operator * ( double s );
    //! Product: Vector b = Matrix A * Vector v
    Vector operator * ( const Vector& v );
    //! Product: Matrix C = Matrix A * Matrix B
    Matrix operator * ( const Matrix& B );
    //----- division operators
    //! Matrix B = Matrix A / s: element-wise division with a scalar
    Matrix operator / ( double s );
    //! Vector b = Matrix A_inv * Vector v: product of inverse matrix and a vector
    Vector operator / ( const Vector& v );
    //! Matrix C = Matrix A_inv * Matrix B: product of inverse matrix and a matrix
    Matrix operator / ( const Matrix& B );
    // inversion and transposition
    //! inverse of Matrix SELF
    const Matrix operator ! ();
    //! transpose of Matrix SELF
    const Matrix operator ~ ();
    //!------------------- console output
    //! prints entries on a single line
    void Print(void);
    //! prints entries line by line
    void pPrint(void);
    //------------------- file I/O
    //! read from predefined file "MATRIX.txt"
    Matrix Read(void);
    //! write to predefined file  "MATRIX.txt"
    void Write(void);
    //! read form a file
    Matrix ReadFrom( const char* file );
    //! write to a file (L entries in a line)
    void WriteTo( const char* file, int L = 20 );
    //!------------------- stream I/O
    friend ostream& operator << ( ostream& s, const Matrix& M );
    //! classify a matrix into diagonal, tridiagonal, symmetric or general type
    void Classify();

  private:
    // addition and subtraction
    //! Matrix B = Matrix self + Matrix A
    Matrix Add( const Matrix& A );
    //! Matrix B = Matrix self - Matrix A
    Matrix Sub( const Matrix& A );
    // various products
    //! Vector b = Matrix self * Vector v
    Vector Dot( const Vector& v );
    //! Matrix B = Matrix self * double s
    Matrix Dot( double s );
    //! Matrix B = Matrix self * Matrix A
    Matrix Dot( const Matrix& A );
    //----- division operators
    //! Matrix B = Matrix A / s: element-wise division with a scalar
    Matrix Div( double s );
    //! Vector b = Matrix A_inv * Vector v: product of inverse matrix and a vector
    Vector Div( const Vector& v );
    //! Matrix C = Matrix A_inv * Matrix B: product of inverse matrix and a matrix
    Matrix Div( const Matrix& B );
    // inversion and transposition
    //! inverse of Matrix SELF
    void Inv();
    //! transpose of Matrix SELF
    const Matrix Tra();
    //! true if matrix is diagonal
    bool isDiag();
    //! true if matrix is tridiagonal
    int isTriDiag();
    //! true if matrix is symmetrix, but not diagonal or tridiagonal
    bool isSymm();
};

//=============================== EXTERNAL MATRIX FUNCTIONS

//! Multiplication of Matrix and a double scalar
Matrix Dot( double s, Matrix& M );
//! Multiplication of Matrix and a double scalar
Matrix Dot( Matrix& M, double s );
//! Multiplication of two matrices: C = A * B
Matrix Dot( Matrix& A, Matrix& B );
//! inverse of Matrix A
Matrix Inv( Matrix& A );
//! transpose of Matrix A
Matrix Tra( Matrix& A );
//! determinant of Matrix A
double Det( Matrix& A );
//! output stream
std::ostream& operator << ( ostream& s, const Matrix& M );
//! LU decomposition for tridiagonal matrices (symmetric or non-symmetric)
void lu_triDiag( Matrix& A );

}  // end of namespace smv_lib

#endif   // _ARRAYS_H

//================================ end of file
